Script to make a first check of the EVA ReSurvey database for MOTIVATE

Load libraries

library(here)
library(readr)
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(stringr)
library(purrr)
library(readxl)

Read the data

db_resurv<-read_tsv(here("data", "edited", "db_resurv.csv"))

Aviso: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)Rows: 425310 Columns: 69── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (31): Country, Biblioreference, Nr. table in publ., Cover abundance scal...
dbl (22): PlotObservationID, PlotID, TV2 relevé number, Nr. relevé in table,...
lgl (16): Lon1, Lon2, Lon3, Lon4, Lat1, Lat2, Lat3, Lat4, X1, X2, X3, X4, Y1...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Problems (do not affect us)

problems<-problems(db_resurv)
sort(unique(problems$col))

[1]  7 13

names(db_resurv[c(7,13)])

[1] "Nr. relevé in table" "Altitude"           

No problems really!

Update coordinates

Create new column with old coordinates if new not available, and with new if available.

db_resurv <- db_resurv %>%
  mutate(Lon_updated = ifelse(is.na(Lon_prec),Longitude,Lon_prec),
         Lat_updated = ifelse(is.na(Lat_prec),Latitude,Lat_prec))

print(db_resurv, width = Inf)

ISSUE 1: ReSurvey plot is NA

Careful! Sometimes ReSurvey plot is NA.

nrow(db_resurv%>%filter(is.na(RS_CODE)))

[1] 0

nrow(db_resurv%>%filter(is.na(`ReSurvey site`)))

[1] 0

nrow(db_resurv%>%filter(is.na(`ReSurvey plot`)))

[1] 111

write_csv(db_resurv %>% filter(is.na(`ReSurvey plot`)),
          here("output", "csv","issue1.csv"))

ISSUE 2: Coordinates are NA

Coordinates are also NA in some cases.

nrow(db_resurv%>%filter(is.na(Lon_updated) & is.na(Lat_updated )))

[1] 949

write_csv(db_resurv %>% filter(is.na(Lon_updated) & is.na(Lat_updated)),
          here("output", "csv","issue2.csv"))

But Resurvey plot and coordinates are never missing at the same time.

nrow(db_resurv%>%
       filter(is.na(`ReSurvey plot`) & is.na(Lon_updated) & is.na(Lat_updated)))

[1] 0

ISSUE 3: Different ReSurvey observations within the same plot have different coordinates

Careful! In some cases different ReSurvey observations within the same plot have different coordinates.

Create two new columns in db_resurv: coordinates_equal indicating if coordinates are exactly equal between ReSurvey observations, and coordinates_consistent, indicating if coordinates are consistent between ReSurvey observations (consistent meaning that difference < 0.001 degrees).

# Define a threshold (e.g., 0.001 degrees for longitude/latitude differences)
threshold <- 0.001

db_resurv <- db_resurv %>%
  group_by(RS_CODE, `ReSurvey site`, `ReSurvey plot`) %>%
  mutate(
    lon_range = ifelse(all(is.na(Lon_updated)), NA,
                        max(Lon_updated, na.rm = T) - 
                         min(Lon_updated, na.rm = T)),
    lat_range = ifelse(all(is.na(Lat_updated)), NA,
                        max(Lat_updated, na.rm = T) - 
                         min(Lat_updated, na.rm = T)),
    coordinates_equal = ifelse(is.na(Lon_updated) & is.na(Lat_updated), NA,
                               lon_range == 0 & lat_range == 0),
    coordinates_consistent = ifelse(is.na(Lon_updated) & is.na(Lat_updated), NA,
                                    lon_range < threshold & 
                                      lat_range < threshold)
  ) %>%
  ungroup() %>%
  select(-lon_range, -lat_range)

write_csv(db_resurv %>% filter(coordinates_equal==FALSE),
          here("output", "csv","issue3.csv"))

db_resurv %>% 
  group_by(RS_CODE,`ReSurvey site`, `ReSurvey plot`) %>%
  summarize(is_equal = all(coordinates_equal),
            is_consistent = all(coordinates_consistent),
            .groups = "drop") %>%
  mutate(coordinate_status = case_when(
    is_equal ~ "Equal",
    !is_equal & is_consistent ~ "Consistent (< 0.001º)",
    !is_equal & !is_consistent ~ "Inconsistent (> 0.001º)")) %>%
  count(coordinate_status)%>%
  mutate(percentage = n / sum(n) * 100) %>%
  ggplot(aes(x = percentage, y = coordinate_status, fill = coordinate_status)) +
  geom_bar(stat = "identity") + 
  geom_text(aes(label = paste0(round(percentage, 1), "%")),
            position = position_stack(vjust = 0.5), size = 3) + 
  labs(x = "Percentage of Plots", y = NULL) +
  theme(axis.text.y = element_text(size = 12)) +
  coord_flip() + theme(legend.position = "none")
ggsave(filename=here("output", "figures","issue3.tiff"),
       width=10,height=7,units="cm",dpi=300)

ISSUE 4: Some plots have only one resurvey

When ReSurvey plot is not NA, use the unique combination of RS_CODE, ReSurvey site and ReSurvey plot to uniquely define each ReSurvey plot (as defined in metadata). When ReSurvey plot is NA (111 rows), use the unique combination of RS_CODE, ReSurvey site and updated coordinates to uniquely define each ReSurvey plot. Check how many resurveys (i.e. different years are there for each unique combination).

count_resurveys <- db_resurv %>%
  # Convert dates to date format and get the year
  mutate(date = dmy(`Date of recording`), year = year(date)) %>%
  group_by(RS_CODE, `ReSurvey site`,
           # If ReSurvey plot is not NA, 
           # group by RS_CODE, `ReSurvey site`, `ReSurvey plot`
           `ReSurvey plot` = ifelse(is.na(`ReSurvey plot`), 
                                    NA_character_, `ReSurvey plot`),
           # If ReSurvey plot is NA, group by coordinates
           Lon_updated = ifelse(is.na(`ReSurvey plot`), Lon_updated, NA_real_),
           Lat_updated = ifelse(is.na(`ReSurvey plot`) , Lat_updated, NA_real_)
  ) %>%
  summarise(
    # Get how many different years for each unique group
    distinct_years=n_distinct(year), 
    # Get how many different dates for each unique group
    distinct_dates=n_distinct(date), .groups = "drop")

Summary stats:

summary(count_resurveys$distinct_years)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   2.000   3.449   4.000  55.000 

sd(count_resurveys$distinct_years)

[1] 2.677649

Histograms:

# For all data
ggplot(count_resurveys, aes(x = distinct_years)) + 
  geom_histogram(fill = "white", color = "black", bins = 55)+
  xlab("Number of ReSurvey observations (different years)") +
  ylab("Number of plots")
ggsave(filename=here("output", "figures","issue4.tiff"),
       width=11,height=7,units="cm",dpi=300)

Number and proportion of plots with only 1 resurvey (should not be so!)

nrow(count_resurveys%>%filter(distinct_years==1))

[1] 1158

nrow(count_resurveys%>%filter(distinct_years==1))/nrow(count_resurveys)

[1] 0.009573571

write_csv(count_resurveys%>%filter(distinct_years==1),
          here("output", "csv","issue4.csv"))

ISSUE 5: Datasets with only presence/absence

db_resurv %>%
  filter(`Cover abundance scale`=="Presence/Absence") %>%
  distinct(Dataset)

ggplot(db_resurv %>% 
         mutate(pres_or_ab =ifelse(`Cover abundance scale`=="Presence/Absence",
                                   "Presence/Absence", "Abundance"),
                DK_Naturdata_Res = ifelse(Dataset == "DK_Naturdata_Res",
                                          "Y", "N")),
                aes(pres_or_ab, fill = DK_Naturdata_Res)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage", x = NULL)
ggsave(filename=here("output", "figures","issue5.tiff"),
       width=12,height=7,units="cm",dpi=300)

For DK_Naturdata_Res - info about habitat from Jerker’s file (see below).

ISSUE 6: Observations with wrong country (GIS)

Read text file wrong_countries obtained in ArcGIS:

wrong_countries <- read_delim(here("data", "clean","wrong_countries.txt"),
                              delim = ";")

Rows: 1850 Columns: 19── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (12): Country, RS_CODE, ReSurvey_s, ReSurvey_p, Lon_update, Lat_update, ...
dbl  (7): FID, Join_Count, TARGET_FID, plot_uniqu, year, obs_unique, PlotObsID
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

write_csv(wrong_countries, here("output", "csv","issue6.csv"))

ISSUE 7: Different cover abundance scales

ggplot(db_resurv, aes(`Cover abundance scale`)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations", x = "Cover abundance scale") +
  coord_flip()
ggsave(filename=here("output", "figures","issue7.tiff"),
       width=18,height=10,units="cm",dpi=300)

ISSUE 8: Wrong EUNIS codes

Used this info in metadata file:

Expert system classification to EUNIS habitats (https://zenodo.org/records/4812736 ; https://floraveg.eu/habitat/). I am sending you legend for EUNIS classification version 2022-10-16 with all codes and meanings, directly prepared from expert system file (second sheet) - it is slightly different from published version in ZENODO (https://zenodo.org/records/4812736 , little bit old dated now) and from https://floraveg.eu/habitat/ (little bit newer than in current EVA version).

Qa = mires and Qb = wetlands P units – in floraveg.eu there is slightly different classification (https://floraveg.eu/habitat/overview/P), but in EVA is still this classification of P:

P Surface waters Pa Base-poor spring and spring brook Pb Calcareous spring and spring brook Pc Brackish-water vegetation Pd Fresh-water small pleustophyte vegetation Pe Fresh-water large pleustophyte vegetation Pf Fresh-water submerged vegetation Pg Fresh-water nymphaeid vegetation Ph Oligotrophic-water vegetation Pi Dystrophic-water vegetation Pj Stonewort vegetation

Presence of “!” simply means that for one unit there are two or more different formulas, e.g. R11 and R11!. So it is only technical stuff.

Multiple assignment of relevé – no priority, alphabetical order, e.g. N16!,S66,S81 means that relevé can be assigned to all 3 units: N16 Mediterranean and Macaronesian coastal dune grassland (grey dune), S66 Mediterranean halo-nitrophilous scrub and S81 Canarian xerophytic scrub

No value present in Expert System – relevé didn´t enter expert system classification (= it means that some prerequisites are missing)

“~” – relevé entered expert classification however was not classified to any EUNIS unit +

Clean info on Expert system column and separate it when there are several codes.

db_resurv <- db_resurv %>%
  mutate(
    # Clean 'Expert System' column by removing "!" and replacing "~" with NA
    `Expert System` = case_when(
      `Expert System` == "~" ~ NA_character_,  # Replace "~" with NA
      TRUE ~ str_replace_all(`Expert System`, "!", "")  # Remove "!"
    )
  ) %>%
  # Separate the values in 'Expert System' into multiple columns
  separate(
    `Expert System`,
    into = c("EUNISa", "EUNISb", "EUNISc", "EUNISd"),
    sep = ",",
    extra = "drop",  # Drop extra values if there are more than columns
    fill = "right",   # Fill missing values with NA for cases with fewer values
    remove = FALSE    # Keep the original 'Expert System' column
  )

Calculate how many different EUNIS codes have been assigned:

db_resurv <- db_resurv %>%
  mutate(
    # Count the number of non-NA values across the EUNIS columns
    n_EUNIS = rowSums(!is.na(select(., starts_with("EUNIS"))))
  )

ggplot(db_resurv, aes(n_EUNIS)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Number of differnt EUNIS codes assigned") + coord_flip()

ggplot(db_resurv %>% filter(n_EUNIS > 0), aes(n_EUNIS)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Number of differnt EUNIS codes assigned") + coord_flip()

Correct some EUNIS codes that are probably wrong:

db_resurv <- db_resurv %>%
  mutate(across(starts_with("EUNIS"), ~ case_when(
    . == "N16M" ~ "N16",
    . == "Sa" ~ "V4",
    . == "Sb" ~ "V5",
    . == "T1CT" ~ "T1C",
    . == "N15A" ~ "N15",
    TRUE ~ .
  )))

Add columns for the different EUNIS levels:

db_resurv <- db_resurv %>%
  mutate(
    # EUNISa levels
    EUNISa_1 = substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 2, 1)),
    EUNISa_2 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 3, 2), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 3, 2)),
      NA_character_
    ),
    EUNISa_3 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 4, 3), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 4, 3)),
      NA_character_
      ),
    EUNISa_4 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 5, 4), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 5, 4)),
      NA_character_
    ),
    
    # EUNISb levels
    EUNISb_1 = substr(EUNISb, 1, ifelse(str_starts(EUNISb, "MA"), 2, 1)),
    EUNISb_2 = ifelse(
      nchar(EUNISb) >= ifelse(str_starts(EUNISb, "MA"), 3, 2), 
      substr(EUNISb, 1, ifelse(str_starts(EUNISb, "MA"), 3, 2)),
      NA_character_
    ),
    EUNISb_3 = ifelse(
      nchar(EUNISb) >= ifelse(str_starts(EUNISb, "MA"), 4, 3), 
      substr(EUNISb, 1, ifelse(str_starts(EUNISb, "MA"), 4, 3)),
      NA_character_
    ),
    EUNISb_4 = ifelse(
      nchar(EUNISb) >= ifelse(str_starts(EUNISb, "MA"), 5, 4), 
      substr(EUNISb, 1, ifelse(str_starts(EUNISb, "MA"), 5, 4)),
      NA_character_
    ),
    
    # EUNISc levels
    EUNISc_1 = substr(EUNISc, 1, ifelse(str_starts(EUNISc, "MA"), 2, 1)),
    EUNISc_2 = ifelse(
      nchar(EUNISc) >= ifelse(str_starts(EUNISc, "MA"), 3, 2), 
      substr(EUNISc, 1, ifelse(str_starts(EUNISc, "MA"), 3, 2)),
      NA_character_
    ),
    EUNISc_3 = ifelse(
      nchar(EUNISc) >= ifelse(str_starts(EUNISc, "MA"), 4, 3), 
      substr(EUNISc, 1, ifelse(str_starts(EUNISc, "MA"), 4, 3)),
      NA_character_
    ),
    EUNISc_4 = ifelse(
      nchar(EUNISc) >= ifelse(str_starts(EUNISc, "MA"), 5, 4), 
      substr(EUNISc, 1, ifelse(str_starts(EUNISc, "MA"), 5, 4)),
      NA_character_
    ),
    
    # EUNISd levels
    EUNISd_1 = substr(EUNISd, 1, ifelse(str_starts(EUNISc, "MA"), 2, 1)),
    EUNISd_2 = ifelse(
      nchar(EUNISd) >= ifelse(str_starts(EUNISd, "MA"), 3, 2), 
      substr(EUNISd, 1, ifelse(str_starts(EUNISd, "MA"), 3, 2)),
      NA_character_
    ),
    EUNISd_3 = ifelse(
      nchar(EUNISd) >= ifelse(str_starts(EUNISd, "MA"), 4, 3), 
      substr(EUNISd, 1, ifelse(str_starts(EUNISd, "MA"), 4, 3)),
      NA_character_
    ),
    EUNISd_4 = ifelse(
      nchar(EUNISd) >= ifelse(str_starts(EUNISd, "MA"), 5, 4), 
      substr(EUNISd, 1, ifelse(str_starts(EUNISd, "MA"), 5, 4)),
      NA_character_
    )
  )

Create new columns with descriptions for the level 1 codes:

db_resurv <- db_resurv %>%
  mutate(
    EUNISa_1_descr = case_when(
      EUNISa_1 == "V" ~ "Vegetated man-made habitats",
      EUNISa_1 == "U" ~ "Inland habitats with no or little soil",
      EUNISa_1 == "T" ~ "Forests and other wooded land",
      EUNISa_1 == "S" ~ "Heathlands, scrub and tundra",
      EUNISa_1 == "R" ~ "Grasslands",
      EUNISa_1 == "Q" ~ "Wetlands",
      EUNISa_1 == "P" ~ "Inland waters",
      EUNISa_1 == "N" ~ "Coastal habitats",
      EUNISa_1 == "MA" ~ "Marine habitats",
      TRUE ~ NA_character_
    ),
    EUNISb_1_descr = case_when(
      EUNISb_1 == "V" ~ "Vegetated man-made habitats",
      EUNISb_1 == "U" ~ "Inland habitats with no or little soil",
      EUNISb_1 == "T" ~ "Forests and other wooded land",
      EUNISb_1 == "S" ~ "Heathlands, scrub and tundra",
      EUNISb_1 == "R" ~ "Grasslands",
      EUNISb_1 == "Q" ~ "Wetlands",
      EUNISb_1 == "P" ~ "Inland waters",
      EUNISb_1 == "N" ~ "Coastal habitats",
      EUNISb_1 == "MA" ~ "Marine habitats",
      TRUE ~ NA_character_
    ),
    EUNISc_1_descr = case_when(
      EUNISc_1 == "V" ~ "Vegetated man-made habitats",
      EUNISc_1 == "U" ~ "Inland habitats with no or little soil",
      EUNISc_1 == "T" ~ "Forests and other wooded land",
      EUNISc_1 == "S" ~ "Heathlands, scrub and tundra",
      EUNISc_1 == "R" ~ "Grasslands",
      EUNISc_1 == "Q" ~ "Wetlands",
      EUNISc_1 == "P" ~ "Inland waters",
      EUNISc_1 == "N" ~ "Coastal habitats",
      EUNISc_1 == "MA" ~ "Marine habitats",
      TRUE ~ NA_character_
    ),
    EUNISd_1_descr = case_when(
      EUNISd_1 == "V" ~ "Vegetated man-made habitats",
      EUNISd_1 == "U" ~ "Inland habitats with no or little soil",
      EUNISd_1 == "T" ~ "Forests and other wooded land",
      EUNISd_1 == "S" ~ "Heathlands, scrub and tundra",
      EUNISd_1 == "R" ~ "Grasslands",
      EUNISd_1 == "Q" ~ "Wetlands",
      EUNISd_1 == "P" ~ "Inland waters",
      EUNISd_1 == "N" ~ "Coastal habitats",
      EUNISd_1 == "MA" ~ "Marine habitats",
      TRUE ~ NA_character_
    )
  )

Plot for EUNISa_1 (the first assigned EUNIS in cases of multiple assignations, level 1):

ggplot(db_resurv, aes(EUNISa_1_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 1") + coord_flip()

ggplot(db_resurv %>% filter(!is.na(EUNISa_1_descr)), aes(EUNISa_1_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 1") + coord_flip()
ggsave(filename=here("output", "figures","issue8.tiff"),
       width=18,height=10,units="cm",dpi=300)

ISSUE 9: Manipulated plots and info on manipulation type

ggplot(db_resurv, aes(`Manipulate (y/n)`)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Manipulation")
ggsave(filename=here("output", "figures","issue9.tiff"),
       width=10,height=8,units="cm",dpi=300)

List of Type of Manipulation in manipulated plots (mixed information):

write_csv(data.frame(unique(db_resurv$`Type of manipulation`)),
          here("output", "csv","issue9.csv"))

ISSUE 10: Location method

ggplot(db_resurv, aes(`Location method`)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Location method") + coord_flip()
ggsave(filename=here("output", "figures","issue10.tiff"),
       width=18,height=8,units="cm",dpi=300)

ISSUE 11: Resurvey project types

unique(db_resurv$RS_PROJTYP)

[1] "resampling"      "permanent"       "permanent (man)" "Permanent (man)"
[5] "Resampling"     

Unify codes:

db_resurv <- db_resurv %>%
  mutate(RS_PROJTYP = recode(RS_PROJTYP,
                             "Resampling" = "resampling",
                             "Permanent (man)" = "permanent (man)"))

unique(db_resurv$RS_PROJTYP)

[1] "resampling"      "permanent"       "permanent (man)"

ggplot(db_resurv, aes(RS_PROJTYP, fill=`Manipulate (y/n)`)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Resurvey project type") + coord_flip() +
  theme(legend.position = "top")
ggsave(filename=here("output", "figures","issue11.tiff"),
       width=18,height=8,units="cm",dpi=300)

ISSUE 12: Column RS_DUPL

db_resurv %>% filter(!is.na(RS_DUPL)) %>% select(RS_CODE, RS_DUPL) %>%
  distinct()

ISSUE 13: Location uncertainty

db_resurv <- db_resurv %>%
  # Redefine precision_new, which was wrong
  mutate(precision_new = factor(ifelse(is.na(Lon_prec) & is.na(Lat_prec),
                                       0, 1)))

ggplot(db_resurv, aes(`Location uncertainty (m)`, fill = precision_new)) +
  geom_histogram( color = "black") +
  xlab("Location uncertainty (m)")

ggplot(db_resurv %>% filter(`Location uncertainty (m)` <= 500),
       aes(`Location uncertainty (m)`, fill = precision_new)) +
  geom_histogram(color = "black") +
  xlab("Location uncertainty (m) <= 500")
ggsave(filename=here("output", "figures","issue13_1.tiff"),
       width=18,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(`Location uncertainty (m)` > 500),
       aes(`Location uncertainty (m)`, fill = precision_new)) +
  geom_histogram(color = "black") +
  xlab("Location uncertainty (m) > 500")
ggsave(filename=here("output", "figures","issue13_2.tiff"),
       width=18,height=8,units="cm",dpi=300)

NO ISSUES FROM HERE

Altitude and slope values

Unique slope values:

unique((db_resurv)$`Slope (°)`) %>% str_sort()

 [1] "_"  "-"  "-1" "."  "0"  "0." "00" "03" "07" "1"  "1." "10" "11" "12" "13"
[16] "14" "15" "16" "17" "18" "19" "2"  "2." "20" "21" "22" "23" "24" "25" "26"
[31] "27" "28" "29" "3"  "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "4" 
[46] "4." "40" "41" "42" "43" "44" "45" "46" "47" "48" "5"  "5." "50" "51" "52"
[61] "55" "58" "6"  "6." "60" "65" "7"  "70" "75" "77" "78" "8"  "8." "80" "85"
[76] "9"  "9." "90" "95" NA  

Set altitude, slope and aspect as numeric:

db_resurv <- db_resurv %>%
  mutate(
    # Some altitude values have a "-" after the number,
    # convert to numeric after removing that
    Altitude = as.numeric(gsub("-", "", Altitude)),
    # Some slope values are noted as "_" or "-", these should be NA,
    # otherwise convert to numeric
    `Slope (°)` = ifelse(`Slope (°)` == "_" | `Slope (°)` == "-",
                   NA, as.numeric(`Slope (°)`)),
    # Convert aspect values to numeric
    `Aspect (°)` = as.numeric(`Aspect (°)`)
    )

Aviso: There was 1 warning in `mutate()`.
ℹ In argument: `Slope (°) = ifelse(`Slope (°)` == "_" | `Slope (°)` == "-", NA,
  as.numeric(`Slope (°)`))`.
Caused by warning in `ifelse()`:
! NAs introducidos por coerción

Histograms:

ggplot(db_resurv, aes(Altitude)) +
  geom_histogram(fill = "white", color = "black")

ggplot(db_resurv, aes(`Aspect (°)`)) +
  geom_histogram(fill = "white", color = "black")

ggplot(db_resurv, aes(`Slope (°)`)) +
  geom_histogram(fill = "white", color = "black")
ggsave(filename=here("output", "figures","issue8.tiff"),
       width=18,height=10,units="cm",dpi=300)

range(db_resurv$`Slope (°)`, na.rm=T)

[1] -1 95

Add columns date and year

db_resurv <- db_resurv %>%
  mutate(date = dmy(`Date of recording`), year = year(date))

Histograms:

ggplot(db_resurv, aes(year)) + geom_histogram(fill = "white", color = "black")

Plot size

ggplot(db_resurv, aes(`Relevé area (m²)`)) +
  geom_histogram(fill = "white", color = "black")

Observations with no info on plot size:

nrow(db_resurv %>% filter(is.na(`Relevé area (m²)`)))

[1] 19599

Cover values (total, trees, shrubs, herbs, mossess)

db_resurv %>%
  pivot_longer(cols = c(`Cover total (%)`, `Cover tree layer (%)`,
                        `Cover shrub layer (%)`, `Cover herb layer (%)`,
                        `Cover moss layer (%)`),
               names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = value)) +
  geom_histogram(fill = "white", color = "black", bins = 10) +
  facet_wrap(~ variable, scales = "free") +
  labs(x = "Value", y = "Frequency")

db_resurv %>%
  reframe(across(c(`Cover total (%)`, `Cover tree layer (%)`,
                     `Cover shrub layer (%)`, `Cover herb layer (%)`,
                     `Cover moss layer (%)`), ~range(., na.rm = TRUE)))

All values OK.

Mosses and lichens identified

ggplot(db_resurv, aes(`Mosses identified (Y/N)`)) + geom_bar()

ggplot(db_resurv, aes(`Lichens identified (Y/N)`)) + geom_bar()

NA in most cases.

All resurveys for each resurvey plot to send to Bea

db_Europa<- db_resurv %>%
  group_by(RS_CODE, `ReSurvey site`,
           # If ReSurvey plot is not NA, 
           # group by RS_CODE, `ReSurvey site`, `ReSurvey plot`
           `ReSurvey plot` = ifelse(is.na(`ReSurvey plot`), 
                                    NA_character_, `ReSurvey plot`),
           # If ReSurvey plot is NA, group by coordinates
           # Create a unique grouping variable that uses coordinates
           # only when conditions are met
           group_coords = ifelse(is.na(`ReSurvey plot`),
                                 paste(Lon_updated, Lat_updated), NA_character_)
  ) %>%
  # Add unique identifiers for each plot.
  # These are based on the unique combination of RS_CODE, ReSurvey site and 
  # ReSurvey plot (When ReSurvey plot is not NA)
  # and on the unique combination of RS_CODE, ReSurvey site 
  # and updated coordinates (When ReSurvey plot is NA)
  mutate(plot_unique_id = cur_group_id()) %>%
  select(PlotObservationID, Country, `Date of recording`, RS_CODE,
         `ReSurvey site`, `ReSurvey plot`, Lon_updated, Lat_updated,
         group_coords, `Expert System`, `Location method`, plot_unique_id) %>%
  ungroup() %>%
  # Convert dates to date format and get the year
  mutate(date = dmy(`Date of recording`), year = year(date)) %>%
  select(-`Date of recording`, -date, -group_coords) %>%
  # Add unique identifiers for each observation
  mutate(obs_unique_id = row_number())

print(db_Europa, width = Inf)

Save to csv (file for us):

write_csv(db_Europa,here("data", "clean","db_Europa_20250107.csv"))

Save to csv (file for Bea, with only essential info):

write_csv(db_Europa %>% 
            select(obs_unique_id, plot_unique_id, Lon_updated, Lat_updated,
                   year),
          here("data", "clean","db_Europa_20241210_short.csv"))

Info on HabitatID from DK

Based on information got from Jesper.

Read the data sent by Jesper from DK

db_DK_J<-read_tsv(here("data", "raw",
                       "DK_Naturdata_Res_habitat_hab_codes_Jesper",
                  "DK_Naturdata_Res_habitat_hab_codes.txt"))

Rows: 158800 Columns: 9── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): HABITAT, Dataset
dbl (7): PlotObservationID, RELEVE_NR, CIRC_ID, PLOT5_ID, PLOT15_ID, Access ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Add info on HabitatID to db_resurv

db_resurv <- db_resurv %>%
  # Keeping all obs in db_resurv but not all in db_DK_J
  left_join(db_DK_J %>% select(PlotObservationID, HabitatID))

Joining with `by = join_by(PlotObservationID)`

List of HabitatID

print(db_resurv %>% distinct(HabitatID), n = 100)

Write csv:

write_csv(db_resurv %>% distinct(HabitatID),
          here("data", "clean","list_HabitatID_DK_resurv.csv"))

Cases without HabitatID info

Cases without ESy EUNIS habitat:

nrow(db_resurv %>% filter(is.na(`Expert System`)))/nrow(db_resurv)

[1] 0.6814982

Cases without ESy EUNIS habitat but with HabitatID from DK:

nrow(db_resurv %>% filter(is.na(`Expert System`)&!is.na(HabitatID)))/nrow(db_resurv)

[1] 0.2727117

Cases without ESy EUNIS habitat and without HABITAT from DK:

nrow(db_resurv %>% 
       filter(is.na(`Expert System`)&is.na(HabitatID)))/nrow(db_resurv)

[1] 0.4087865

Cases without ESy EUNIS habitat and without HabitatID from DK where data is presence / absence:

nrow(db_resurv %>%
       filter(is.na(`Expert System`) &
                is.na(HabitatID) &
                `Cover abundance scale` == "Presence/Absence"))/
  nrow(db_resurv)

[1] 0.1655851

Cases without ESy EUNIS habitat and without HabitatID from DK where data is not presence / absence:

nrow(db_resurv %>%
       filter(is.na(`Expert System`) &
                is.na(HabitatID) &
                `Cover abundance scale` != "Presence/Absence"))/
  nrow(db_resurv)

[1] 0.2432014

Change some Annex I habitat codes that were wrong

db_resurv <- db_resurv %>%
  mutate(HabitatID = as.character(HabitatID)) %>%
  mutate(HabitatID = ifelse(HabitatID == "9998", "91D0",
                            ifelse(HabitatID == "9999", "91E0", HabitatID)))

Add info on correspondences HabitatID (DK, Jesper) - EUNIS

Read correspondences file:

correspondences<-read_excel(here("data", "edited",
                                 "correspondence_HabitatID_DK.xlsx"))

Add info to db_resurv:

db_resurv <- db_resurv %>%
  # Keeping all obs in db_resurv but not all in db_DK_J
  left_join(correspondences %>% select(HabitatID, EUNIS))

Joining with `by = join_by(HabitatID)`

Correct NA values in EUNIS

db_resurv <- db_resurv %>%
  mutate(EUNIS = ifelse(EUNIS == "NA", NA, EUNIS))

Add info on EUNIS (DK) to EUNISa:

db_resurv <- db_resurv %>%
  mutate(EUNISa =
           # If EUNIS (DK) is available, add as EUNISa
           ifelse(!is.na(EUNIS), EUNIS, 
                  # Otherwise keep EUNISa
                  EUNISa),
         EUNIS_assignation = ifelse(!is.na(EUNIS), "Info from DK",
                                    ifelse(is.na(EUNISa), "Not possible",
                                           "Expert system"))) %>%
  # Remove column EUNIS (DK)
  select(-EUNIS)

ggplot(db_resurv, aes(EUNIS_assignation)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS assignation")

Generate file for ISSUE 5

write_csv(db_resurv %>%
            filter(is.na(`Expert System`)&is.na(HabitatID)),
          here("output", "csv","issue5.csv"))

Update columns for EUNIS levels and descriptions

Update the columns for the different EUNISs levels:

db_resurv <- db_resurv %>%
  mutate(
    # EUNISa levels
    EUNISa_1 = substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 2, 1)),
    EUNISa_2 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 3, 2), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 3, 2)),
      NA_character_
    ),
    EUNISa_3 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 4, 3), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 4, 3)),
      NA_character_
      ),
    EUNISa_4 = ifelse(
      nchar(EUNISa) >= ifelse(str_starts(EUNISa, "MA"), 5, 4), 
      substr(EUNISa, 1, ifelse(str_starts(EUNISa, "MA"), 5, 4)),
      NA_character_
    )
  ) %>%
  # Remove HabitatID column
  select(-HabitatID)

Update columns with descriptions for the level 1 codes:

db_resurv <- db_resurv %>%
  mutate(
    EUNISa_1_descr = case_when(
      EUNISa_1 == "V" ~ "Vegetated man-made habitats",
      EUNISa_1 == "U" ~ "Inland habitats with no or little soil",
      EUNISa_1 == "T" ~ "Forests and other wooded land",
      EUNISa_1 == "S" ~ "Heathlands, scrub and tundra",
      EUNISa_1 == "R" ~ "Grasslands",
      EUNISa_1 == "Q" ~ "Wetlands",
      EUNISa_1 == "P" ~ "Inland waters",
      EUNISa_1 == "N" ~ "Coastal habitats",
      EUNISa_1 == "MA" ~ "Marine habitats",
      TRUE ~ NA_character_
    )
  )

Number of different EUNIS codes

Recalculate how many different EUNIS codes have been assigned:

db_resurv <- db_resurv %>%
  mutate(
    # Count the number of non-NA values across the EUNIS columns
    n_EUNIS = rowSums(!is.na(select(., EUNISa:EUNISd)))
  )

ggplot(db_resurv, aes(n_EUNIS)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Number of differnt EUNIS codes assigned") + coord_flip()

ggplot(db_resurv %>% filter(n_EUNIS > 0), aes(n_EUNIS)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "Number of differnt EUNIS codes assigned") + coord_flip()

New plot for EUNISa_1 (the first assigned EUNIS in cases of multiple assignations, level 1):

ggplot(db_resurv, aes(EUNISa_1_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 1") + coord_flip()

ggplot(db_resurv %>% filter(!is.na(EUNISa_1_descr)), aes(EUNISa_1_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 1") + coord_flip()

Add info on descriptions for EUNIS levels 2-4

descriptions<-read_excel(here("data", "edited",
                                 "EUNIS-Habitats-2021-06-01_modified.xlsx"))

# Define the columns and corresponding description column names
eunis_cols <- c("EUNISa_2", "EUNISa_3", "EUNISa_4",
                "EUNISb_2", "EUNISb_3", "EUNISb_4", 
                "EUNISc_2", "EUNISc_3", "EUNISc_4",
                "EUNISd_2", "EUNISd_3", "EUNISd_4")

# Create corresponding description column names
descr_col_names <- paste0(eunis_cols, "_descr")

# Use reduce to loop through the columns and join dynamically based on level
db_resurv <- reduce(seq_along(eunis_cols), function(data, i) {
  # Extract level number from the column name (e.g., EUNISa_2 -> 2)
  level <- as.numeric(gsub("\\D", "", eunis_cols[i]))
  
  # Filter descriptions for the corresponding level
  descriptions_level <- descriptions %>%
    filter(level == level) %>%
    select(`EUNIS 2020 code`, `EUNIS-2021 habitat name`)
  
  # Perform the left_join and rename the column dynamically
  data %>%
    left_join(
      descriptions_level,
      by = setNames("EUNIS 2020 code", eunis_cols[i])
    ) %>%
    rename(!!descr_col_names[i] := `EUNIS-2021 habitat name`)
}, .init = db_resurv)

The matching did not work sometimes, correct!

# Correct EUNISa levels 2-4 descriptions
db_resurv <- db_resurv %>%
  mutate(EUNISa_2_descr = 
           ifelse(!is.na(EUNISa_2_descr), EUNISa_2_descr,
                  case_when(
                    EUNISa_2 == "Pf" ~ "Fresh-water submerged vegetation",
                    EUNISa_2 == "Pj" ~ "Stonewort vegetation",
                    EUNISa_2 == "R4" ~ "Alpine and subalpine grasslands",
                    EUNISa_2 == "Pb" ~ "Calcareous spring and spring brook",
                    EUNISa_2 == "Qb" ~ "Wetlands",
                    EUNISa_2 == "R3" ~ "Seasonally wet and wet grasslands",
                    EUNISa_2 == "Qa" ~ "Mires",
                    EUNISa_2 == "Pa" ~ "Base-poor spring and spring brook",
                    EUNISa_2 == "Ph" ~ "Oligotrophic-water vegetation",
                    EUNISa_2 == "Pg" ~ "Fresh-water nymphaeid vegetation",
                    EUNISa_2 ==
                      "Pd" ~ "Fresh-water small pleustophyte vegetation",
                    EUNISa_2 == "Pc" ~ "Brackish-water vegetation",
                    EUNISa_2 ==
                      "Pe" ~ "Fresh-water large pleustophyte vegetation",
                    EUNISa_2 == "Pi" ~ "Dystrophic-water vegetation",
                    EUNISa_2 == "S1" ~ "Tundra",
                    EUNISa_2 ==
                      "U7" ~ "Unvegetated or sparsely vegetated gravel bars",
                    EUNISa_2 == "Q6" ~ "Periodically exposed shores",
                    TRUE ~ NA_character_)
                  ),
         EUNISa_3_descr = 
           ifelse(!is.na(EUNISa_3_descr), EUNISa_3_descr,
                  case_when(
                    EUNISa_3 =="U71" ~ "Unvegetated or sparsely vegetated gravel bar in montane and alpine regions",
                    EUNISa_3 =="Q61" ~ "Periodically exposed shore with stable, eutrophic sediments with pioneer or ephemeral vegetation",
                    EUNISa_3 =="Q62" ~ "Periodically exposed shore with stable, mesotrophic sediments with pioneer or ephemeral vegetation",
                    TRUE ~ NA_character_
                    ))
         )

# Correct EUNISb levels 2-4 descriptions
db_resurv <- db_resurv %>%
  mutate(EUNISb_2_descr = 
           ifelse(!is.na(EUNISb_2_descr), EUNISb_2_descr,
                  case_when(
                    EUNISb_2 == "Pj" ~ "Stonewort vegetation",
                    EUNISb_2 == "R4" ~ "Alpine and subalpine grasslands",
                    EUNISb_2 == "Pf" ~ "Fresh-water submerged vegetation",
                    TRUE ~ NA_character_)
                  )
         )

EUNISc and EUNISd levels 2-4 are OK.

Notes EUNIS codes - to change?

https://www.sci.muni.cz/botany/chytry/Schaminee_etal2021_EEA-Report-Aquatic-Wetland-habitats.pdf

EUNISa_2 == “Q6” : “Periodically exposed shores” EUNISa_3 = “Q61” : “Periodically exposed shore with stable, eutrophic sediments with pioneer or ephemeral vegetation” EUNISa_3 == “Q62” : “Periodically exposed shore with stable, mesotrophic sediments with pioneer or ephemeral vegetation”

This classification of Q + numbers is now coexisting in the database with Qa & Qb (metadata). How to proceed?

db_resurv %>% filter(EUNISa_1 == "Q") %>% distinct(EUNISa_2)

Plots of level-2 categories within each level 1 category

ggplot(db_resurv %>% filter(EUNISa_1 == "MA"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "MA") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","MA_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "N"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "N") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","N_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "P"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "P") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","P_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "Q"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "Q") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","Q_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "R"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "R") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","R_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "S"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "S") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","S_level2.tiff"),
       width=16,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "T"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "T") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","T_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "U"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "U") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","U_level2.tiff"),
       width=16,height=8,units="cm",dpi=300)

ggplot(db_resurv %>% filter(EUNISa_1 == "V"), aes(EUNISa_2_descr)) +
         geom_bar(aes(y = (..count..) / sum(..count..) * 100)) +
  labs(y = "Percentage of ReSurvey observations",
       x = "EUNIS level 2") + coord_flip() +
  ggtitle(db_resurv %>% filter(EUNISa_1 == "V") %>% distinct(EUNISa_1_descr))
ggsave(filename=here("output", "figures","V_level2.tiff"),
       width=14,height=8,units="cm",dpi=300)

Save to clean data

Save so-far clean datafile for resurvey database:

write_tsv(db_resurv,here("data", "clean","db_resurv_clean.csv"))

Session info

sessionInfo()

R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
[3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
[5] LC_TIME=Spanish_Spain.utf8    

time zone: Europe/Paris
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] readxl_1.4.3    purrr_1.0.2     stringr_1.5.1   ggplot2_3.5.1  
[5] lubridate_1.9.4 dplyr_1.1.4     tidyr_1.3.1     readr_2.1.5    
[9] here_1.0.1     

loaded via a namespace (and not attached):
 [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    stringi_1.8.4    
 [5] hms_1.1.3         digest_0.6.37     magrittr_2.0.3    evaluate_1.0.1   
 [9] grid_4.4.2        timechange_0.3.0  fastmap_1.2.0     cellranger_1.1.0 
[13] rprojroot_2.0.4   jsonlite_1.8.9    fansi_1.0.6       scales_1.3.0     
[17] jquerylib_0.1.4   cli_3.6.3         rlang_1.1.4       crayon_1.5.3     
[21] bit64_4.5.2       munsell_0.5.1     withr_3.0.2       cachem_1.1.0     
[25] yaml_2.3.10       tools_4.4.2       parallel_4.4.2    tzdb_0.4.0       
[29] colorspace_2.1-1  vctrs_0.6.5       R6_2.5.1          lifecycle_1.0.4  
[33] bit_4.5.0.1       vroom_1.6.5       pkgconfig_2.0.3   pillar_1.9.0     
[37] bslib_0.8.0       gtable_0.3.6      glue_1.8.0        xfun_0.49        
[41] tibble_3.2.1      tidyselect_1.2.1  knitr_1.49        farver_2.1.2     
[45] htmltools_0.5.8.1 rmarkdown_2.29    labeling_0.4.3    compiler_4.4.2